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Abstract 

We adapt the Kolmogorov's normalization algorithm (which is the key element 
of the original proof scheme of the KAM theorem) to the construction of a suitable 
normal form related to an invariant elliptic torus. As a byproduct, our procedure can 
also provide some analytic expansions of the motions on elliptic tori. By extensively 
using algebraic manipulations on a computer, we explicitly apply our method to a 
planar four-body model not too different with respect to the real Sun-Jupiter- 
Saturn-Uranus system. The frequency analysis method allows us to check that our 
location of the initial conditions on an invariant elliptic torus is really accurate. 
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1 Introduction 

Since the birth of the KAM theory (see [20], [36] and [TJ), the invariant tori are expected 
to be the key dynamical object which explains the (nearly perfect) quasi-periodicity of 
the planetary motions of our Solar System. 

Among the consequences of the KAM theory, which concern the tori of maximal 
dimension, the following one looks natural. In a planetary system including n planets, 
one expects that the persistence under small perturbations should hold also for the n- 
dimensional invariant tori, which are a slight deformation of the composition of n coplanar 
circular Keplerian orbits. However, a separate proof is needed in order to ensure the 
existence of these lower dimensional invariant tori which are said to be elliptic, because 
they correspond to stable equilibrium points of the secular motions. Such a theorem has 
been recently proved by Biasco, Chierchia and Valdinoci in two different cases: for the 
spatial three-body planetary problem and for a planar system with a central star and 
n planets (see [I] and [5], respectively). In our opinion, their approach is deep from 
a theoretical point of view, but is not suitable for explicit applications, even if one is 
interested just in finding the locations of the elliptic invariant tori. In order to clarify 
this point, let us roughly summarize the scheme of their proofs as follows: first, they 
carry out all the preliminary canonical transformations that are necessary to bring the 
Hamiltonian in a particular form, to which they can subsequently apply a theorem due 
to Poschel (see [ID] and [41]). so to ensure the existence of elliptic lower dimensional tori. 
Moreover, Poschel's versions of this theorem are based on a careful adaptation of the 
usual Arnold's proof scheme for non-degenerate systems: the perturbation is removed by 
a sequence of canonical transformations which are defined on a subset of the phase space 
excluding the "resonant regions" (see [T] and [2]). Since resonances are everywhere dense 
(but the width of the regions eliminated around them is suitably decreased, when the 
order of the resonances increases), therefore the change of coordinates giving the shape 
of the invariant elliptic tori is defined on a Cantor set which does not contain any open 
subset. The efficiency of an eventual explicit application based of such an approach is 
highly questionable and, as far as we know, it has never been used to calculate an orbit 
of a Celestial Mechanics problem. 

The original proof scheme of the KAM theorem, introduced by Kolmogorov himself, 
is in a much better position for what concerns the translation into an explicit algorithm 
constructing invariant tori (see [20], [3], [13] and JH]). In fact, this approach has been suc- 
cessfully used to calculate the orbits for some interesting problems in Celestial Mechanics 
(see [32], [33], [31] and [H]). The present work aims to adapt the Kolmogorov's algorithm, 
in order to construct a suitable normal form related to the elliptic tori. Moreover, this 
will allow us to explicitly integrate the equations of motion on those invariant surfaces, 
by using a so called semi-analytic procedure. 

When one is interested in showing the long term stability of a planetary system, the 
construction of a normal form related to some fixed elliptic torus could be a relevant 
milestone. In fact, it is possible to ensure the effective stability in the neighborhood of 
such an invariant surface by implementing a partial construction of the Birkhoff 's normal 
form (see, e.g., [18] and [15], where this approach is used in order to study the stability 
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nearby an invariant KAM torus having maximal dimension). For what concerns our 
Solar System, such an approach might be applied to some asteroids with small orbital 
eccentricities and inclinations. However, as explained in [13], this same approach can not 
yet succeed in proving the long-time stability of the major planets of our Solar System. 

The location of the elliptic tori can be useful also for practical purposes. In fact, 
the regions close to them are exceptionally stable, being mainly filled by invariant tori 
of maximal dimension. Therefore, they can be of interest for spatial missions aiming, 
for instance, to observe asteroids not far from the elliptic tori. Moreover, our technique 
should adapt quite easily also to the construction of hyperbolic tori that can be used in 
the design of spacecraft missions with transfers requiring low energy. Also in view of this 
kind of applications, lower dimensional tori of elliptic, hyperbolic and mixed type have 
been studied in the vicinity of the Lagrangian points for both the restricted three-body 
problem and the bicircular restricted four-body problem (see, e.g., [17], [19], [8] and [TO]). 

The present paper is organized as follows. The search for elliptic tori is applied just 
to a model not far from the SJSU planar system (let us recall that the real orbits of the 
planets of our Solar System are not lying on lower dimensional tori). Therefore, sect. |2] 
is devoted to the introduction of our Hamiltonian model and to the description of its 
expansion in canonical coordinates. This will allow us to write down the form of the 
Hamiltonian to which our approach can be applied. By the way, we think that with some 
minor modifications our procedure should adapt also to the more general spatial case, after 
having performed the reduction of the angular momentum, which is not considered here in 
order to shorten the description of all the preliminary expansions (for an introduction to 
some methods performing both the partial and the total reduction, see [9], [35] and [37]). 

Our algorithm constructing a normal form for elliptic tori is presented in a purely 
formal way in sect. |3j Let us recall that our procedure is mainly a reformulation of the 
classical Kolmogorov's normalization algorithm, that is modified in a suitable way for our 
purposes. The theoretical background necessary to understand when our algorithm can 
converge is informally discussed in subsect. 13.41 and it will be fully detailed in a future 
work (see [T6]). 

Sect. |H is devoted to explain an application which is also a test of our procedure. 
First, in subsects. I4.1H4.2I we describe the way to implement our algorithm, by using an 
algebraic manipulator on a computer so to produce both the normal form and the semi- 
analytic integration of the motion on an invariant elliptic torus. The Fourier spectrum of 
the motions on elliptic tori is strongly characteristic: just the mean- motion frequencies 
and their linear combinations can show up. This simple remark allows us to check the 
accuracy of our results by using frequency analysis, as it will be described in sect. 14.31 

2 Classical expansion of the planar planetary Hamilto- 
nian 

As claimed in the introduction, in order to fix the ideas, we think it is convenient to focus 
on a concrete planetary model, to which we will apply our algorithm constructing the 
elliptic tori in the next sections. Let us consider four point bodies Pq, Pi, P2, P3, with 
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masses m , mi, m 2 , m 3 , mutually interacting according to Newton's gravitational law. 
For shortness, hereafter we will assume that the indexes 0, 1, 2, 3 correspond^ to Sun, 
Jupiter, Saturn and Uranus, respectively. 

Let us now recall how the classical Poincare variables can be introduced so to perform 
a first expansion of the Hamiltonian around circular orbits, i.e., having zero eccentricity. 
We basically follow the formalism introduced by Poincare (see [3B] and [21]; for a modern 
exposition, see, e.g., [27] and [28]) . We remove the motion of the center of mass by 

using heliocentric coordinates r ■ =PqPj , with j = 1, 2, 3 . Denoting by r_- the momenta 
conjugated to r,-, the Hamiltonian of the system has 6 degress of freedom, and reads 



F(r , r) = T (0) (£) + U {0) (r) + T (1) (f) + U (1) (r) , (1) 



where 



3 
I 
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The plane set of Poincare's canonical variables is introduced as 



T<°>(r) = lE^lfell 2 , r«(r) = 4(r 1 -f 2 + r 1 .f 3 + f 2 .f 3 ) 

j=i \ 

a i 3 " \ "— 1 — 2 " "—1 —3 II 



3 > 



mg W3 



m m 7 - , , , , 



(2) 



for j = 1,2,3, where cij , , Mj and cjj are the semi-major axis, the eccentricity, 
the mean anomaly and the perihelion argument, respectively, of the j—th planet. One 
immediately sees that both £j and rjj are of the same order of magnitude as the eccentricity 
Cj . Using Poincare's variables (E}, the Hamiltonian F can be rearranged so that one has 

F(A,\,i,v) = ^ (0) (A) + F {1 \A,X,i,v) , (3) 

where = T^+U^, F« = T^ + U^. Let us emphasize that = 0(1) and F« = 
0(/i) , where the small dimensionless parameter /x = max{m! / m , m 2 / m , m 3 / m } 
highlights the different size of the terms appearing in the Hamiltonian. Therefore, let us 
remark that the time derivative of each coordinate is O(fi) but in the case of the angles 
A . Thus, according to the common language in Celestial Mechanics, in the following we 
will refer to A and to their conjugate actions A as the fast variables, while (£, rj) will be 
called secular variables. 

We proceed now by expanding the Hamiltonian ([3]) in order to construct the first basic 
approximation of the normal form for elliptic tori. After having chosen a center value A* 



-"^Let us stress that the four considered point bodies have the same masses as Sun, Jupiter, Saturn and 
Uranus, but the orbits here studied are significantly different with respect to the real ones. 
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Table 1: Masses rrij and initial conditions for Jupiter, Saturn and Uranus in our planar 
model. We adopt the UA as unit of length, the year as time unit and set the gravitational 
constant Q = 1 . With these units, the solar mass is equal to (27r) 2 . The initial conditions 
are expressed by the usual heliocentric planar orbital elements: the semi-major axis cij , 
the mean anomaly Mj , the eccentricity Cj and the perihelion longitude Uj . The data are 
taken by JPL at the Julian Date 2440400.5 . 





Jupiter (J — 1) 


Saturn (J = 2) 


Uranus (j = 3) 


m j 

CLj 
Mj 

e 3 


(2tt) 2 /1047.355 

5.20463727204700266 

3.04525729444853654 

0.04785365972484999 

0.24927354029554571 


(2tt) 2 /3498.5 

9.54108529142232165 

5.32199311882584869 

0.05460848595674678 

1.61225062288036902 


(2vr) 2 /22902.98 

19.2231635458410572 

0.19431922829271914 

0.04858667407651962 

2.99374344439246487 



for the Taylor expansions with respect to the fast actions (in a way we will explain later), 
we perform a translation 7a* defined as 

L 3 =A 3 -A*, Vj = l,2,3. (4) 

This is a canonical transformation that leaves the coordinates A, £ and r\ unchanged. 
The transformed Hamiltonian HP^ = F o 7a* can be expanded in power series of L, £, rj 
around the origin. Thus, forgetting an unessential constant we rearrange the Hamiltonian 
of the system as 

oo oo oo 

(L, \,Z_,r 1 ) = n*.L + J2 (^) + EE h n)n & A, £, v) , (5) 

ji=2 ji=0j 2 =0 

(T) 

where the functions ,- 2 are homogeneous polynomials of degree ji in the actions L 
and of degree j'2 in the secular variables (£, rj) . The coefficients of such homogeneous 

polynomials do depend analytically and periodically on the angles A. The terms h^ e ^ 
of the Keplerian part are homogeneous polynomials of degree ji in the actions L, the 
explicit expression of which can be determined in a straightforward manner. In the latter 
equation the term which is both linear in the actions and independent of all the other 
canonical variables (i.e., n* ■ L) has been separated in view of its relevance in perturbation 
theory, as it will be discussed in the next section. We also expand the coefficients of the 
power series h^ F j 2 in Fourier series of the angles A. The expansion of the Hamiltonian 
is a traditional procedure in Celestial Mechanics. We work out these expansions for the 
case of the planar SJSU system using a specially devised algebraic manipulator. The 
calculation is based on the approach described in sect. 2.1 of [32J, which in turn uses the 
scheme sketched in sect. 3.3 of [42J. 

The reduction to the planar case is performed as follows. We pick from Table IV of jB] 
the initial conditions of the planets in terms of heliocentric positions and velocities at the 
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Julian Date 2440400.5 . Next, we calculate the corresponding orbital elements with respect 
to the invariant plane (that is perpendicular to the total angular momentum). Finally we 
include the longitudes of the nodes Qj (which are meaningless in the planar case) in the 
corresponding perihelion longitude Uj and we eliminate the inclinations by setting them 
equal to zero. The remaining initial values of the orbital elements are reported in Table [U 

Having fixed the initial conditions we come to determining the average values (a^ , a 2 , a%) 
of the semi-major axes during the evolution. To this end we perform a long-term numer- 
ical integration of Newton's equations starting from the initial conditions related to the 
data reported in Table [TJ After having computed (a\ , a 2 , a%) , we determine the val- 
ues A* via the first equation in ([2]). This allows us to perform the expansion of the 
Hamiltonian as a function of the canonical coordinates (£, A, £,77). In our calculations 
we truncate this initial expansion as follows, (a) The keplerian part is expanded up to 
the quartic terms. The series where the general summand - 2 appears are truncated 
so to include: (bl) the terms having degree j\ in the actions L with ji < 3 , (b2) all 
terms having degree j 2 in the secular variables (£, rj) , with j 2 such that 2j 1 + j 2 < 8 , 
(b3) all terms up to the trigonometric degree 18 with respect to the angles A. Let us 
remark that with respect to the analogous initial expansion we performed in [43], here 
we preferred to considerably reduce the maximal degree in the secular coordinates, in 
order to increase those related to the fast ones. This choice is motivated by the fact that 
the orbits on elliptic tori experience smaller values of the eccentricities (let us recall that 
both £j = 0(ej) and rjj = 0(ej) V j = 1 , 2 , 3 ) than those related to the real motions; 
moreover, larger limits on the fast coordinates are needed, in order to give a sharp enough 
numerical evidence of the convergence of the algorithm described in the next section. 

Let us now focus on the average with respect to the fast angles of the Hamiltonian 
written in §5§, i.e. ("H ( - r ^) A . The fast actions L are obviously invariant with respect 

to the flow of (T-L ( - T ') x , thus, they can be neglected while just the secular motions are 
considered. The remaining most significant term is given by the lowest order approxim- 
ation of the secular Hamiltonian, namely its quadratic term (/102 ) \ > which is essentially 
the one considered in the theory first developed by Lagrange (see j2TJ) and furtherly im- 
proved by Laplace (see [21], [23] and [26]) and by Lagrange himself (see [22], [23J). In 
modern language, we can say that the origin of the secular coordinates phase space (i.e., 
= (0)0) ) is an elliptic equilibrium point for the secular Hamiltonian. In fact, under 
mild assumptions on the quadratic part of the Hamiltonian which are satisfied in our case 
(see sect. 3 of [5], where such hypotheses are shown to be generically fulfilled for a planar 
model of our Solar System), it is well known that one can find a canonical transformation 
(L, A, £,77) = V(p,q,x,y) owning the following properties: (i) L — p and A = q, (ii) the 
map (£,77) = ((,(x),v(y)) is linear, (iii) V diagonalizes the quadratic part of the Hamilto- 
nian, so that we can write (h^ 2 ) x in the new coordinates as Y^j=i u j (rf + Vj)/^ > where 
all the entries of the vector u} ^ have the same sign. Our algorithm constructing a suitable 
normal form for elliptic tori can be started from the Hamiltonian = HP* o T> , i.e. 

H^(p,q,x,y) = H^(V(i,q,x,y)) . (6) 
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3 Formal algorithm 

In the present section, let us more generically assume that the number of degrees of 
freedom of our system is n\ + U2 , where the canonical coordinates (p, q, x, y) can naturally 
be split in two parts, that are (p, q) E W 11 x TP 1 and (x,y) E W 12 xl" 2 ." 

In order to better understand our whole procedure, we think it is convenient to 
immediately state our final goal. We want to determine a canonical transformation 
(P,q,x,y) = /C (oo) (P,Q,X,F} such that the Hamiltonian #(°°) = #(°) o /C (oo) is brought 
to the following normal forno 

ri2 q(oo) / y2 | -ia2^ 

(P, Q, X, Y) = -E+Yl - [ J j) + 

3=1 2 ( ?) 

o(iipii 2 ) + o(iipiiii(x,i:)||) + o(ii(x,i:)|| 3 ) , 

where the notation means that we want to remove all terms which are linear in P_ and 
independent of ( X , Y_) , or at most quadratic in (X, Y_) and independent of P. 

When initial conditions of the type (P,Q,X,Y) = (0,Q ,0,0) (with Q Q E T ni ) are 
considered, the normal form (J2J) allows us to easily calculate the solution of the Hamilton 
equations, i.e. 

(P(t),Q(t),X(t),Y(t)) = (0,g o + u; (oo) t,0,0) . (8) 

This clearly means that the rii-dimensional (elliptic) torus corresponding to P = X = 
Y_ = is invariant and the orbits are quasi-periodic on it. 

Let us start the description of the generic r-th step of our algorithm constructing the 
normal form. We begin with a Hamiltonian of the following type: 

oo oo 

H (r - 1] (g, q, x, y) = <Jr-*> ■ p + ■ J + E fci? (& ^ & ' ( 9 ) 

s=0 1=0 2j 1 +j 2 =l 

31>0,J 2 > 

where Jj = (x'j + y])/2 is the action which is usually related to the j-th pair of secular 
canonical coordinates (xj,yj) , V j = 1 , . . . , n 2 . Moreover, there is a fixed integer value 
K > such that the terms ff^j^ 1 satisfy the following hypotheses: 

(A) fjljz E Vjlfj 2 , where Vjlj 2 is the class of functions such that (al) they are 
homogeneous polynomials of degree j± in the actions p , (a2) they are homogeneous 
polynomials of degree j'2 m the secular variables (x, y) , (a3) their Fourier expansion 
is finite with maximal trigonometric degree equal to sK ; 

(B) the terms fjl~j 2 ' s ^ are "well Fourier-ordered" ; this nonstandard definition means that 
V ji > , J2 > , s > 1 every Fourier harmonic k appearing in the expansion of 



2 Let us here stress a little abuse of notation. Hereafter, the symbol ui will mean the frequencies vector 
related to the motion on a torus (as it is usual in KAM theory), while in the previous sections it was 
used to represent the perihelion longitudes (according to the classical notation in Celestial Mechanics). 
Analogously, hereafter, will denote the oscillation frequencies transverse to an elliptic torus, while 
before it was used for the longitudes of the nodes. 
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•mi h * s sucn ^ na t ^ s corresponding trigonometric degree \k\ = |fcx| + . . . + \kj^\ > 
(s- l)K . 

By using formula (j5J) and the properties (i)-(iii) of the canonical transformation T> , 
one easily sees that the Hamiltonian (that is defined in (jSJ)) can be expanded in the 
form written in after having suitably reordered its Fourier expansion so to satisfy the 
above requirements (A) and (B). Therefore, our constructive algorithm can be concretely 
applied to the Hamiltonian by starting with r = 1 . 

The comparison of the expansion in with the normal form in ([7j) clearly shows 
that we have to eliminate all the terms fj^j 2 where the index I = 2ji + is such that 
< I < 2 . Thus, the r-th step of our algorithm can be naturally divided in three stages, 
each of ones aims to reduce the perturbative terms with I = 0, 1, 2, respectively. 

3.1 First stage of the r-th normalization step: removing of the 
terms depending just on q 

By making use of the classical Lie series algorithm to calculate canonical transforma- 
tions (see, e.g., [12] for an introduction), we first introduce the new Hamiltonian H^ I;r ' = 
exp£ (rjif ( r-1 \ where the generating function Xo (q) ^00 i s determined as the solu- 
tion of the equation 

{xP,^.p} + j2tio 1 ' S \l) = 0, (10) 

8=1 

where we used the classical symbol {-, •} to represent the Poisson brackets. The previous 
equation (that is usually said to be of homological type) admits a solution provided the 
frequency vector w^ -1 -* is non-resonant up to order rK , i.e. 

min \k ■ a/ 1 " -1 - 1 1 > a r with a r > , (11) 

0<|fc|<r_K' 1 1 

where, for the time being, {a r } r> o is nothing but a sequence of real positive numbers and 
\k\ denotes the / 1 -norm of the integer vector k, i.e. \k\ = \ki \ + . . . + |fe m | . The solution 
of the homological equation (JTUJ) can be easily recovered by looking at the little more 
complicate case of , which is discussed in the third stage of the r-th normalization 
step (see formulas (j2T|) - (j23]) ). 

In order to avoid the proliferation of too many symbols, let us make a common abuse of 
notation so to still denote with (p, q, x, y) the new canonical coordinates exp C^ r ) (p, q, x, y) . 
The expansion of the new Hamiltonian can be written as follows: 

00 00 

ircw (g f gj £ y j = jt-i) . £ + a (r-D E fit* fe> £• a y) ■ ( 12 ) 

S=0 1=0 2h+32=t 
J1>0,J2>0 

The mathematical recursive definitions of the terms ff^ r f 2 are lenghty, but it is rather 
easy to understand how to deal with them when they are translated in a programming 
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language. The main remark is concerned with the classes of functions, i.e. 

U'vffc^ G Pf+SP V0<K Jllj2 >0, S >0. (13) 

Therefore, after having calculated all the Poisson brackets needed by the term |,£ l ( r )fj[~j^ S ' 

it is enough to know that it contributes to the sum X^=o fji-i,h ' ^ suitable "reordering 
of the Fourier series" will allow us to ensure that also the expansion ()12p satisfies the 
conditions (A) and (B), which have been stated at the beginning of the present section. 



3.2 Second stage of the r-th normalization step: removing of 
the terms linear in (x, y) and independent of p 

Let us now introduce the new Hamiltonian if ( II;r ) = exp C ( T )H^ 1I ' r \ where the generating 

(r) / \ (rK) 

function x\ {<liX,y) G Vq A is determined as the solution of the equation 



ii2 n( r_1 ) 



(14) 



s=0 



In order to explicitly write down the solution of the previous equation, it is convenient to 
temporarily introduce action-angle coordinates so to replace the secular pairs (x,y) by 
putting Xj = ^/2Jj cos (fj and jjj = sin cpj V j — 1, . . . , n-i ; therefore, let us assume 

that the expansion of the known terms appearing in equation ()14p is the following one: 



0<\k\<rK j=l 



s=0 



with suitable real coefficients c^j and . Thus, one can easily check that 



II 2 



0<|fc|<rK 3=1 



j. sin (fc • g ± J cos (k_ ■ q ± ipj) 



(15) 



(16) 



is a solution of the homological equation (j!4p and it exists provided the frequency vector 
u/ r ~ 1 ) satisfies the so-called first Melnikov non-resonance condition up to order rK , i.e. 



mm 

0<|fe|<r_ftT 

3=1 , ... ,712 



with a r > , 



(17) 



and all the entries of the frequency vector ^ are far enough from the origin, i.e. 



min \Q 

j=l,...,n 2 



(r-l)| 



> /3 with > 



For what concerns planetary Hamiltonians where the D'Alembert rules hold true, let us 

,3 H,3 



remark that all the coefficients cj^ ] and d^- appearing in (ITS"]) and having even values of 
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\k\ are equal to zero. In order to solve the equation ( THI) . therefore, we could to not need 
the condition fll8p . which however is substantially included in another one (i.e., (I3T|) ) that 
we will be forced to introduce later. 

Starting from the expansion ffTBT) of Xi (9S J-i V 9 ) > one can immediately recover the 
expression of Xi (<7>:& u) as a function of the original polynomial variables. We can then 
explicitly calculate the expansion of the new Hamiltonian, which can be written as follows: 

oo oo 

H^ r \p,q,x,y)=^- 1) -p + Q^-J + Y,Y, E fj^h&l'&y) ■ (19) 

S=0 1=0 2J1+J 2 =' 
31>0,J 2 >0 

Also in this case, providing mathematical recursive definitions of the terms is a 

quite ennoying task. Thus, we think it is better to just describe how to deal with them 
when they are translated in a programming language. First, let us remark that the 
following relations about the classes of functions hold true: 

T&p E ftS e U v nT ] V < i < I , s > . (20) 

2ji+j2=l 2j 1 +j 2 =l-i 

Therefore, after having calculated all the Poisson brackets appearing in the expression 
of the term j^C\ r) ^2 2 ji+j 2 =i fjl'^ja '> ^ * s enou gh to know that it contributes to the sum 

5^=o Yl2j 1 +j 2 =i-i fju?^ • Again a suitable "reordering of the Taylor-Fourier series" will 
allow us to ensure that also the expansion fll9p satisfies the conditions (A) and (B), which 
have been stated at the beginning of the present section. 



3.3 Third stage of the r-th normalization step: removing both 
the linear terms in p that are independent of (x, y) and the 
quadratic ones in (x, y) which are independent of p 

The Hamiltonian produced at the end of the r-th normalization step is provided by the 
composition of three canonical transformation^] which can be given in terms of Lie series, 
i.e. = exp£ (r) o exp£ (r) o exp C ( r )H^ u ' r \ where the generating functions belong 

to three different classes: X% (p, q) € , ^ 2 (q,x,y) e ^02 an d (x,y) G Vq \ . 

The explicit expressions of these generating functions are given below, in formulas (T2~Tj) . 
(121) and respectively. 

3 When one focuses on the estimates needed to prove the convergence of the algorithm, it is certainly 
simpler to introduce the generating function x£ (p,Q,x,y) — X^~\p,q) + (q,x,y) and to consider 
the new Hamiltonian exp£_(r> o exp£ (r)-ff < - II;r ) which slightly differs from H^ r \ because X% > and 

(r) 

Y2 do not commute with respect to the Poisson brackets. Since in the present work we do not want to 
theoretically study the problem of the convergence of our algorithm, wc think here is better to distinguish 
the present third stage of the r-th normalization step in other three parts, so to highlight their different 
roles. Moreover, this choice looks more natural to us, when one implements the constructive algorithm 
by algebraic manipulations on a computer. 
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We start with (p, q) 6 v[ r ^ , which is determined as the solution of the equation 
{xi r \^.p} + j2f§ r ' S) (P,l) = U- (21) 

This implies that 



X. 



(r) 



(p,g)= E E^ 

0<\k\<rK j=l 



Ckj sin (k ■ q) d^j cos (k ■ q) 

k ■ a/'" 1 ) + k ■ U^" 1 ) 



(22) 



where we preliminarly assumed that the expansion of the known terms appearing i 
equation (12T|) has the form 



in 



E /£) (&£) = E E^' fe'j cos (- ■ i) + d ^>j sin (- ■ s)l ' ( 23 ) 

0<\k\<rK i=l 



s=l 



with suitable real coefficients Ckj and dkj . Let us here recall that the solution f l22|) for 
the equation (T2"Tj) exists provided the frequency vector u/ r_1 ) satisfies the non-resonance 
condition ([Til . 

Let us now consider K 2 (<7> ^ 2/) e ^02 > which is a solution of the following homolo- 
gical equation: 



m n( r '- 1 ) "I r 

Yf ] , ^ -p-Y^i'l • $ + E /S ir,s) (£, £ y) = 



(24) 



8=1 



In order to explicitly write down the expansion of , it is convenient to temporarily 
reintroduce the action-angle coordinates (J, <p) so to replace the secular pairs (x, y) ; 
therefore, let us assume that the expansion of the known terms appearing in equation 02"" 
has the form 



»2 



Y,fS$' r ' s) k,iS= E E 2 v^ 

s=l 0<|&|<rX i,j = l 



(25) 



with suitable real coefficients c^- ^ and dj^' . Thus, one can easily check that 



o<|fc|<rif «,i = i 



+ 



k • c^-d ± of _1) ± n ( ;- 1] 



(26) 
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is a solution of equation ( I2"lj) and it exists provided the frequency vector lo} t ^ satisfies 
the so-called second Melnikov non-resonance condition up to order rK , i.e. 



mm 

o<\k\<rK 

i, j=l , ... ,n 2 



\k ■ u {r - l) ± _1) ± nS P_1) I > a r with a r > . (27) 



Let us here remark that the previous assumption includes also the non-resonance con- 
dition ffTlj) as a special case, i.e. when i = j and the signs appearing in the expression 



±f2j- r ^ ± f^ r ^ are opposite. 

Also for what concerns the generating function once again it is convenient to 
replace the secular pairs (x, y) with the action-angle coordinates (J, (p) . Let us here 

remark that ff r ~^ ■ J and fo^'^i^y) are the only terms appearing in expansion (Tl9|) . 
quadratic in (x, y) (so they also are 0( J) ) and not depending on p and q . The canonical 

transformation induced by the Lie series exp £_(r) aims to eliminate the part of fo£'°^ 

depending on the secular angles ip . Therefore, the generating function is defined so 
to solve the following equation: 



, n^-j} + f^Uv) - </o ( :?'°\ = o , (28) 

where denotes the average with respect to the angles ip . This implies that 



ri2 

,(r) 



^ r) (z^)= E E 



i ,j = 1 Si ,Sj = ±1 
Si-i+Sj-j ^ 



c iijjSijSj sin (s^i + SjV?i) 

+ ajrfr-v 



djj >Si , Sj cos (s^i + Sjifj) 



(29) 



where we preliminarly assumed that the expansion of the known terms appearing in 
equation (|28|) is the following one: 



/o ( ?' W)= E E 2 v^ 



i,j' = l s»=±l 
s,-=±l 



di,j, Sl ,s 3 sin (sj^i + Sj^j) 



(30) 



with suitable real coefficients Cij jSii Sj and dij jSii Sj . Let us remark that the solution (f22|) 
for the equation (|21|) exists provided the frequency vector f2 ( - r ~ 1 ' ) satisfies the following 
finite non-resonance condition: 

min \l_ ■ | > with > . (31) 
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At this point of the algorithm, it is convenient to slightly modify the frequencies u and 
£2, so to include the terms which are linear with respect to the actions, do not depend 
on the angles and, then, can not be eliminated by our normalization procedure. More 
definitely, we define and f£- r ' , so that 

uA) . p = Jr-X) . £ + f (y,0 )(l) tf) . j = . j + (/^^ . ( 32 ) 

Standard utilities provided by any computer algebra system should allow everyone to 
get the expansions of Y 2 \q, x, y) and T> 2 (x, y) , starting from those of Y 2 (q, J, ip) and 

are written in (1261) and (|29|) . respectively. Thus, we are now able to 
explicitly produce the expansion of the new Hamiltonian, which can be written as follows: 

oo oo 

H (r \p,q,x,y)=u^-p + Q^-J + J2J2 E f£1iJhl>&y) • ( 33 ) 

s=0 1=0 aji+M=i 
?i>o,j 2 >o 

Let us remark that this expansion of H^ r ' has exactly the same form of that written 
for if( r_1 ) in 02]), but we stress that the algorithm is arranged so to make smaller and 
smaller the contribution of the terms fj[j 2 , when the value of r is increased, V s > and 
1=31+32 = 0, 1, 2. 

In this case too, we avoid to write down the lenghty mathematical recursive definitions 
of the terms fj[j 2 ■ Instead, we provide some relations about the classes of functions, 
which are useful to understand how to translate this third stage of the r-th normalization 
step in a programming language. For what concerns the generating function \ the 
following relations about the classes of the functions hold true: 

\ C x^ fn^ e V nT ] V < > > 3i > , h > , s > . (34) 
The relations involving the generating function are a little more complicated: 

E f^ft U V nT ] V,>0,/>0,,>0. (35) 

2ii+j 2 =i 2ji+j 2 =« 

Finally, one can easily remark that each class of function is invariant with respect to a 
Poisson bracket with the generating function T>^\ therefore: 

^' n> /;"r G Vgl V i > , j, > , J2 > , S > . (36) 

By taking into account the relations fl34l) - fl36l) about the classes of functions, the defini- 
tion ( 13 2 p of the new frequencies vectors and by suitably "reordering" the Taylor-Fourier 
series, it is possible to ensure that also the expansion fl3"3"|) satisfies the conditions (A) 
and (B), which have been stated at the beginning of the present section about the equa- 
tion ([9]). Therefore, the whole normalization procedure, that has been here described for 
the r-th step can be iteratively repeated. 
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3.4 Some remarks about the convergence of the normalization 
algorithm 

We devote this section to an informal discussion of the relations between the normal- 
ization procedure for an elliptic torus, which is the subject of the present paper, and 
the Kolmogorov's algorithm for a torus of maximal dimension. Our aim is to bring into 
evidence, on one hand, the differences that make the case of an elliptic lower dimensional 
torus definitely more difficult and, on the other hand, the impact that these differences 
have on the explicit calculation. 

The main hypotheses of the classical KAM theorem are (a) that the perturbation 
should be small enough and (b) that a strong non-resonance condition must be satisfied 
by the frequencies of the unperturbed torus. Both these conditions appear also in the 
proof of existence of elliptic tori, but the condition of non-resonance presents some critical 
peculiarities. 

Concerning the smallness of the perturbation, the main problem remains that the 
analytical estimates are extremely restrictive. Nevertheless, one can obtain realistic res- 
ults by using algebraic manipulation in order to implement a computer-assisted proof, 
as it has been made sevelar times about the classical KAM theorem (see, e.g., [32J). 
A computer-assisted procedure takes advantage of the preliminar application of the al- 
gorithm constructing the normal form (which is explicitly performed for a finite number 
of steps R, as large as possible), because a suitable version of the KAM theorem is fi- 
nally applied to the Hamiltonian H^ R ' having the perturbing terms strongly reduced with 
respect to the initial H^ '. In the case of lower dimensional elliptic tori, by comparing 
the Hamiltonian normal form (J7|) with the expansion (133]) of H^ r \ one easily realizes that 
the initial expression of the perturbation (making part of the Hamiltonian H^°\ written 
in flS])) is given by 

oo 2 

EE E fni- < 37 > 

S=0 1=0 2j 1 +j 2 =i 
il>0,J 2 >0 

Looking at all the preliminary expansions, which have been described in sect. [2] and 
allowed us to introduce the initial Hamiltonian H^°\ one immediately sees that all the 
perturbing terms appearing in fl37|) are proportional to fi . Let us also recall that the 
small parameter \i is equal to the mass ratio between the biggest planet and the central 
star (according to its definition given in the discussion following formula (J3])). In the 
present context, the explicit application of the normalization algorithm mainly requires 
to translate in a programming language the method described in the previous sections. 

From a conceptual point of view, a much more difficult problem is here concerned with 
the conditions on non-resonance for the frequencies. In the case of a torus having maximal 
dimension, one must choose the n frequencies Ux, . . . , u n so as to satisfy a strong non- 
resonance condition. A typical request is that they obey a Diophantine condition, i.e., that 
the sequence {a r } r >i appearing in the inequality (ITT)) must be such that a r > j/(rK) T 
with suitable positive values of the constant 7 and r . This choice must be made at the 
very beginning of the procedure, and the perturbed invariant torus that is found at the 
end has the same frequencies as the unperturbed one. The reason is that at every step a 
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small translation of the actions is introduced in order to keep the frequencies constant. 

In the case of the elliptic lower dimensional torus one deals instead with two separate 
set of frequencies, namely G W 11 which characterize the orbits on the torus, and 
the transverse frequencies G W 12 that are related to the oscillation of orbits close to 
but not lying on the torus. By the way, this justifies the adjective "transverse" that is 
commonly used. Now, the frequencies oA ^ on the torus can be chosen in an arbitrary 
manner, but the transverse frequencies are functions of being given by the 
Hamiltonian. This is easily understood by considering the case of a periodic orbit, i.e., 
tlx — 1, since in that case the transverse frequencies are related to the eigenvalues of the 
monodromy matrix. 

The striking fact is that, due precisely to the dependence of the transverse frequencies 
on to} the algorithm forces us to change these frequencies at every step. That is, one 
actually deals with infinite sequences u/- r ^ and all required to satisfy at every order 
a non-resonance condition of the form ( |27p . Moreover, both sequences should converge 
to a final set of frequencies = wHfyM) and fi (oo) = fi (oo) (<JV) which must be 

non-resonant (e.g., Diophantine). Thus, we are forced to conclude that, depending on the 
initial choice of u}°\ it may happen that the algorithm stops at some step because the 
frequencies fail to satisfy at least one of the non-resonance conditions f lTTj) . (IT7|) . ffTHj) . (I27|) 
and (|3ip . This is indeed one of the main difficulties in working out the proof of existence 
of an elliptic torus. 

Let us first consider the analytical aspect, in such a way that we can sketch some 
of the ideas that will be exploited in detail in a future more theoretical work, dedicated 
to the same subject studied here. One initially focus on an open ball B C M ni such 
that the Diophantine condition at finite order required for the first step is satisfied by 
every u/°) G B C IR ni and by the corresponding transverse frequencies This can be 
done, because only a finite number of non-resonance relations are considered. Therefore, 
one shows that at every step there exists a subset of frequencies in B which satisfies 
the non-resonance conditions (still at finite but increasing order) required in order to 
perform the next step, together with the corresponding transverse frequencies. This is 
obtained by a procedure which is reminiscent of Arnold's proof scheme of KAM theorem: 
at every step one removes from B a finite number of intersections of B with a small 
strip around a resonant plane in W 11 , assuring that the width of the strip decreases fast 
enough so that the remaining set always takes a non-empty interior part. By the way, 
this is also strongly reminiscent of the process of construction of a Cantor set. The final 
goal is to prove just that one is left with a Cantor set on non-resonant frequencies which 
satisfy the required resonance conditions and has positive Lebesgue measure. Moreover, 
the relative measure with respect to B tends to 1 when the size of the perturbation is 
decreased to zero. This is the idea underlying the proof that will be expanded in [16]. 
We emphasize that the procedure outlined here is strongly inspired by the proof scheme 
of KAM theorem introduced by Arnold, which is quite different from Kolmogorov's one 
(compare J2U] with [TJ). 

Let us now come to the numerical aspect. At first sight the formal algorithm seems 
to require a cumbersome trial and error procedure in order to find the good frequencies: 
when some of the non-resonance conditions fail to be satisfied at a given step one should 
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change the initial frequencies and restart the whole process. Moreover, since the non- 
resonance condition must be satisfied by the final frequencies, which obviously can not 
be calculated, the whole process seems to be unsuitable for a rigorous proof. We explain 
here in which sense the computer-assisted proofs can help to improve the results also in 
this context. We make two remarks. 

The first remark is connected with the use of interval arithmetic while performing 
the actual construction. Following the suggestion of the analytic scheme of proof, we 
look for uniform estimates on a small open ball B , such that V G B we explicitly 
perform R normalization steps, with R as large as possible. Essentially, we may reproduce 
numerically the process of eliminating step by step the unwanted resonant frequencies by 
suitably determining the intervals. Once R steps have been explicitly performed, we may 
apply to the partially normalized Hamiltonian H^ R > a suitable formulation of the KAM 
theorem for elliptic tori. This means that we recover the scheme that we have already 
applied to the case of tori with maximal dimension. That is, we can take advantage of 
the fact that the perturbing terms are much smaller than the corresponding ones for the 
initial Hamiltonian H^; thus, in principle we could ensure that for realistic values of /x 
the relative measure of the invariant tori is so large that the set of those a/ -* for which 
the algorithm can not work (i.e., B\S) is so small that can be neglected when we are 
dealing with a practical application. 

Taking a more practical attitude, we may rely on the fact that the set of good fre- 
quencies, according to the theory, has Lebesgue measure close to one, so that the case 
of frequencies which are resonant at some finite order occurs with very low probability. 
Thus, we just make a choice of the initial frequencies and proceed with the construction, 
checking at every order that the non-resonance conditions that we need at that order are 
fulfilled. We emphasize that the most extended resonant regions are those of low order, 
so that it is not very difficult to check initially that the chosen frequencies will likely be 
good enough. It may happen, of course, that the whole procedure must be restarted with 
different frequencies, but we expect that this will rarely occur. However, since the size of 
the perturbation is expected to decrease geometrically, we may confidently expect that 
the probability of failure will decrease, too. This is confirmed by the actual calculations. 

When R steps have been made, in principle we can apply the theorem to a small 
neighborhood of the calculated frequencies by choosing a suitable initial ball around the 
frequencies approximated at that step. 

4 Elliptic tori for the SJSU system 

We come now to the application of the formal algorithm for the construction of an elliptic 
torus to the planar SJSU system. 

The initial Hamiltonian is written in fl5]), with a suitable rearrangement of terms so 
that it is given the form fl9]) with r = 1 . This requires also a diagonalization of the 
quadratic part in the secular variables, which is performed as described at the end of 
sect. M 

In the present section, we explicitly construct the normal form at a finite order checking 
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Norms of the generating functions 




1e-07 



1e-08 



4 5 6 

Normalization step: r 



Figure 1: Algorithm constructing the normal form related to an elliptic torus for the 
planar SJSU system: plot of the norm of the generating functions as a function of the 
normalization step r ; more precisely, the symbols x , □ , A , Q an d + refer to the norm 



(r) 



Xi , X%, Y2 an d ^2 > respectively, which are defined 



of the generating functions Xo 
during the normalization algorithm, as described in sect. [3j The norm is calculated by 
simply adding up the absolute values of all the coefficients appearing in the expansion of 
each generating function. 
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that the norms of the generating function decrease as predicted by the theory. Then we 
perform a numerical check by comparing the orbit obtained via the normal form with the 
numerically integrated one. 



4.1 Constructing the elliptic torus by using computer algebra 

We applied the algorithm constructing elliptic tori (which has been widely described 
in sect. [3]) to the Hamiltonian (that is defined in ([6]) and has been obtained as 

described in sect. [2]). The parameters have been fixed according to the specific values 
of the planar SJSU system, which are reported in Table [TJ Our software package for 
computer algebra allowed us to explicitly calculate all the expansions (ED of with 
index r ranging between and 9 , so to include: (cl) the terms having degree j\ in the 
actions p with ji < 3 , (c2) all terms having degree j 2 in the variables (x, y) , with j 2 such 
that 2j 1 + j 2 < 8 , (c3) all terms up to the trigonometric degree 18 with respect to the 
angles q. Let us recall that the truncation rules (cl)-(c3) are in agreement with those 
prescribed about the expansion flS]) in sect. [2] at points (bl)-(b3). Let us remark that 
both the truncation rules (cl) and (c2) are preserved by all the canonical transformations 
included in our algorithm. Moreover, we have found that fixing K = 2 is a suitable 
choice to have a rather regular decreasing of the size of the generating functions when the 
normalization step r is increased, as shown in Fig. [TJ Since the maximal trigonometric 
degree of the generating functions Xo ■> X i ■> X 2 and Y% is equal to rK , the choice to 
set K = 2 and the rule (c3) explain why we stopped the algorithm after having ended the 
normalization step with r = 9 . 

The behavior of the norms of the generating functions is reported in Fig. [TJ Let us 
make a few comments. The theoretical estimates predict that the norms should decrease 
geometrically with the order in order to assure the convergence of the normal form. The 
figure shows that this is indeed the behavior in our case. We emphasize that the presence 
of a dangerous resonance would be reflected in a sudden increase of the coefficients; thus, 
the plot gives a practical confirmation that the frequencies are well chosen. 

Let us stress that performing the construction of the normal form up to order r = 9 has 
been very stressing for the computational resources available to us, although the length 
of the calculation is not reflected in a corresponding length of the present subsection. 



4.2 Explicit calculation of the orbits on the elliptic torus 

We now perform a check on the approximation of the elliptic torus. To this end, we 
calculate the orbit on the torus using the analytic expression and we compare it with 
a numerical integration of Hamilton's equations. In this subsection we explain how the 
calculation of the orbit via normal form is performed. 

According to the theory of Lie series, the canonical transformation (p, q, x, y) = 
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jqt) (p( r \ q( r \ xy\ y( r >) inducing the normalization up to the step r is given by 

expC m o exp C (r) o ... o exp£ « o exp£ wo fag) 

Xl Xo ^2 r 2 v ' 

exp£ (i) o exp C m o exp C m (p (r) , g (r) , x {r) , y {r) ) , 

A 2 Xl Xo v — — — ' 

where {p^ r \ q^ r \ x} r \ y^) are meant to be the new coordinates. Thus, the canonical 
transformation (p,q,x,y) = fC(°°>(P, Q, X, Y) brings in the normal form = 
H'°) o K,(°°\ which is written in (J7J), with /C^ 00 -* = lim^oo KM\ Let us introduce a new 
symbol to denote the composition of all the canonical change of coordinates defined in 
sects. [2] and El i.e. 

CM = s o 71* o V o /C (r) , (39) 

where (r, r) = £(A, A, £, rj) is the canonical transformation giving the heliocentric positions 
r and their conjugated momenta r as a function of the Poincare variables. If (£(0),r(0)) 
is an initial condition on an invariant elliptic torus, in principle we might use the following 
calculation scheme to integrate the equation of motion: 

(c {oo) ) _1 

(£(0),r(0)) — ► (£(0) = , Q(0) , X(0) = 0, Y(0) = 0) 

K(~lp > ( 4 °) 

^(oo) 

(£(*) , r(t)) <— (£(*) = , Q(t) = Q(0) + w (oo) t , X(t) = , F(t) = 0) 

where $*(«,). P induces the quasi-periodic flow related to the frequencies vector u/°°). Of 
course, the previous scheme requires an unlimited computing power; from a practical 
point of view, we can just approximate it, by replacing with C^ R \ where R is as large 
as possible. Thus, the integration via normal form actually reduces to a transformation of 
the coordinates of the initial point to the coordinates of the normal form, the calculation 
of the flow at time t in the latter coordinates, which is a trivial matter since the flow is 
exactly quasi-periodic with known frequencies, followed by a transformation back to the 
original coordinates. 

Such an approximated semi-analytic calculation scheme can be directly compared with 
the results provided by a numerical integrator. As it has been shown in [52], [53], [51] 
and [TTJ, this kind of comparisons provide a very stressing test for the accuracy of the 
whole algorithm constructing the normal form. 
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4.3 Validation of the results by using frequency analysis 

The ideal calculation scheme fHUj) highlights that the Fourier spectrum of each component 
of the motion law t i— > (r(t),r(t)) is the very peculiar one 

oo 

y~] cj exp (iQt) , where, V j > , Cj G C and 3 kj G Z ni such that (j = kj • ^ (oc) • 

3=0 

(41) 

In other words, the Fourier spectrum of the planetary motions on elliptic tori is so charac- 
teristic, because all its frequencies are given by linear combinations of the fast frequencies. 
From a strictly mathematical point of view, let us recall that the previous formula for the 
Fourier spectrum can be deduced by the scheme fj^Ol) . because of the analyticity of the so 
called conjugacy function Q h> fO, Q,0, 01 and this will be ensured as a byproduct 
of the theoretical (future) study of the convergence of the constructive algorithm. 

In the present subsection we aim to check the peculiar quasi-periodicity of the motions 
on our approximation of an elliptic torus, by using the frequency map analysis (see, 
e.g., [2J] and [30] for an introduction). We focus on the following initial conditions: 

(C (9) ) _1 (0,0,0,0) ; (42) 

according to the calculation method described in the previous subsect. I4.2[ this should 
be an accurate approximation of a point on an elliptic torus. Therefore, we preliminarly 
integrated the motion of the planar SJSU system over a time interval of 2 24 years, by 
using the symplectic method SBAB 3 (see [31 J) with a time-step of 0.04 years. 

Here we should add a remark concerning the precision. In order to have a signal 
clean enough to be analyzed a particular care about the precision is mandatory. After 
some trials tuning the parameters of the numerical integration, we found that the 80 bits 
floating point numbers provided by the current AMD and INTEL CPUs fits our needs. 
Technically this is obtained by using the long double types of the GNU C compiler 
under a Linux operating system. 

The orbits have been sampled with a time interval of 1 year. The signals related to 
the secular Poincare variables, that are £/(£) +ir]i(t) with / = 1, 2, 3 , have been submitted 
to the frequency analysis method using the so-called Hanning filter. 

In Tabled we report our numerical results about the first 25 summands of the decom- 
position (j4T|) for the Uranus secular signal, i.e. ^(t) + ir] 3 (t) . Let us point out that the 
values of the fast frequencies vector have been preliminarly calculated by looking 
at the main components of the Fourier spectrum of the signals A;(t)exp (iA;(t)) , with 
/ = 1, 2, 3 . Moreover, we stress that the vectors Ay G Z™ 1 listed in the third column are 
determined so to minimize the absolute difference \Q — kj ■ a/ 00 -*! with \kj\ < 20 ; indeed, 
one has to fix some limits on the absolute value of kj , in order to make consistent its 
calculation, and our choice is motivated by the fact that the Fourier decay of the analytic 
conjugacy function C (oo) (0,Q,0,0) is such that the main contributions to the spectrum 
are related to low order harmonics. 

If the initial conditions ( 14"2]) were exactly on an elliptic torus, each value \Q — kj -a/ 00 )) 
reported in the fourth column of Table [2] should be equal to zero. All of them, except for 
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Table 2: Decomposition of the Fourier spectrum of the signal £ 3 (t) + ii] 3 (t) , which is 
related to the Uranus secular motion. The following numerical values have been obtained 
by applying the frequency analysis method. See the text for more details. 
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Figure 2: Frequency analysis of the secular signal related to the secular Jupiter motion: 
£x(t) +ir]i(t) = Y^JLo c j ex P {Kjt) ■ °f the amplitudes \cj\ as a function of the frequen- 
cies Q in Log-Log scale. The symbol x [+, resp.] refer to the signal related to the motion 
starting from the initial conditions (|42p resp.], i.e. the approximation of a point on 

an elliptic torus after having performed 9 [0, resp.] steps of the algorithm constructing 
the corresponding normal form. In both cases, the results for just the first 25 components 
have been reported in the figure above. 
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the case corresponding to j = 16 , are actually small enough to be considered as generated 
by round-off errors. On the other hand, we can definitely say that £i6 ~ —1.1 x 10~ 5 is a 
"secular frequency", because its value is O(fi) . Indeed, let us recall that /i ~ 10~ 3 , but 
the mass ratio for Uranus, i.e. m^/mo ~ 4.4 x 10 -5 , is even smaller. 

Let us say that the occurrence of secular frequencies in the Fourier decomposition of 
the signal should be expected. Indeed, they could be completely avoided only in a very 
ideal situation, namely: (i) all the calculations described in sects. H] and [3] should be carried 
out without performing any truncations on the expansions, (ii) the initial conditions (j42p 
should be replaced with (C (oo) )~ (0,0,0,0) , (hi) no numerical errors should be there. In 
a practical calculation the orbit can not be exactly placed on an elliptic torus, so the 
presence of secular frequencies just means that we are just close to it. Nevertheless, it is 
very remarkable that the amplitude of the first found secular frequency is three orders of 
magnitude smaller than the main component of the spectrum. In our opinion, this is a 
first clear indication that our algorithm is properly working. 

Other components corresponding to secular frequencies are expected to be even smaller 
than that found with j = 16 . In fact, let us recall that the frequency analysis method 
detects the summands Cj exp (iCj^) appearing in ( 14 ip in a nearly decreasing order with 
respect to the amplitude \cj\ (for instance, one can easily see that just one exchange is 
needed in order to rewrite Table |2]in the correct decreasing order); moreover, we calculated 
that the discrepancy \£ 3 (t) + in 3 (t) — Y^j=o c j ex P (Kjt) | is smaller than about ~ 3.2 x 10~ 7 
for all the time values t for which we sampled the signal. Let us emphasize that such an 
upper bound on the maximal discrepancy is just a little larger than the amplitude | Ci6 1 . 

A similar decomposition has been calculated for both the signals £i(£) + W7i(A) ano - 
^(t) (which are related to the secular motions of Jupiter and Saturn, respectively). 

The behavior is very similar to that of TableEJ so we omit the corresponding tables because 
the interesting results are more evident from the figures that we are going to present. 

The most relevant information about such decompositions of the secular motions of the 
three planets is summarized in the plots done with the x symbol appearing in Figs. [2H11 

Those figures contain also a comparison with the results provided by a, say, trivial 
approximation of an orbit on an elliptic torus. In fact, the dots marked with the + symbol 
appearing in Figs. [2HH refer to a frequency analysis which is performed exactly in the same 
way as that corresponding to the x symbol, except the fact that the numerical integration 
of the equations of motion is started from the following initial conditions 

(C (0) ) _1 (0,0,0,0) , (43) 

instead of that reported in formula ( 142]) . Let us remark that (C^) 1 (0,0,0,0) = £ o 
7a* o T> (0, 0, 0, 0) is a sort of trivial approximation of a point on the elliptic torus as 
it is provided by simply avoiding to apply the part of our algorithm constructing the 
normal form, as it is described in sect. El In order to discuss in a more definite way, 
we have already assumed that the secular frequencies are O(fi) , thus, let us separate 
them from the fast ones, when they are smaller than 10~ 3 . By looking at the right side 
of Figs. EHll one can immediately remark that the parts of the spectra related to the 
fast frequencies are nearly indistinguishable when the initial conditions ( 142]) or ( 143]) are 



24 



M. Sansottera, U. Locatelli, A. Giorgilli 



Saturn frequencies 



0.001 



0.0001 



-B 1e-05 



E 
< 



1e-06 



1e-07 







X 






X 


+ 


+ 


X : 


+ 
















XX 






X %x 






x * x 




X 


X 






X $ 






X 









1e-05 



0.0001 



0.001 0.01 
Frequencies 



0.1 



10 



Figure 3: Frequency analysis of the secular signal related to the secular Saturn motion: 
£2^) + = J2°jLo c j ex P OC?'*) ■ Pl°t °f the amplitudes | Cj I cLS £L function of the fre- 

quencies Q in Log-Log scale. The meaning of the symbols x and + is the same as in 
Fig. El 
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Figure 4: Frequency analysis of the secular signal related to the secular Uranus motion: 
£ 3 (i) + ir]s(t) = J2JLo c j ex P (Kjt) ■ Plot °f the amplitudes \cj\ as a function of the fre- 
quencies (j in Log-Log scale. The meaning of the symbols x and + is the same as in 

Fig. ® 
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considered, because the dots marked with the symbols x and + superpose each other 
in a nearly exact way for what concerns all the main components. On the other hand, 
the secular parts of the spectra (that are in the left side of Figs. [2H4J) strongly differ. In 
fact, when the initial conditions (|43j) (that trivially approximate a point on the elliptic 
torus) are considered, three secular frequencies are detected; while just one is found in 
the case of the more accurate initial data (I42p . Moreover, by comparing the amplitudes, 
one can see that the unique secular component detected by both the frequency analyses 
is decreased by at least two orders of magnitude when our algorithm is applied. In our 
opinion, this comparison makes evident the effectiveness of our procedure constructing 
the normal form for an elliptic torus. 
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